Percentage of Regional Richness

Run ’download_data.Rmd` First!

options(na.action = "na.fail") 

library(tidyverse)
library(dplyr)

source("./helper__dredge_functions.R")

model_average_richness <- function(data) {
  model_average(percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + urban + realm + number_of_checklists_scaled + (1|city_name), data, 4)
}

Find city predictors

create_city_dataset <- function(mixed_model) {
  city_coeefs <- ranef(mixed_model)$city_name
  city_coeefs <- cbind(city_name = rownames(city_coeefs), city_coeefs)
  rownames(city_coeefs) <- NULL
  names(city_coeefs) <- c('name', 'response')
  
  left_join(city_coeefs, city_data)
}

Merlin: Difference in richness

merlin_richness_result <- model_average_richness(merlin)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

merlin_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <8 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/3/4/5/6/7/8/9/10/11    18 -24617.98 49272.05  0.00   0.24
1/2/3/4/5/6/7/8/10/11      17 -24619.11 49272.29  0.24   0.21
1/2/4/5/6/7/8/9/10/11      17 -24619.43 49272.93  0.89   0.15
1/2/4/5/6/7/8/10/11        16 -24620.49 49273.04  0.99   0.14
1/2/3/4/5/6/7/8/9/10/11/12 19 -24618.08 49274.25  2.20   0.08
1/2/3/4/5/6/7/8/10/11/12   18 -24619.09 49274.26  2.22   0.08
1/2/4/5/6/7/8/10/11/12     17 -24620.47 49275.02  2.97   0.05
1/2/4/5/6/7/8/9/10/11/12   18 -24619.52 49275.12  3.07   0.05

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 14.09316    0.88391     0.88404  15.942  < 2e-16 ***
closed_forest                3.89893    0.70680     0.70690   5.516  < 2e-16 ***
cultivated                  11.37400    0.98130     0.98144  11.589  < 2e-16 ***
elevation_delta_scaled       0.58915    0.82679     0.82687   0.713   0.4761    
herbaceous_vegetation        9.88097    1.08890     1.08906   9.073  < 2e-16 ***
herbaceous_wetland          12.35171    3.13811     3.13856   3.935 8.30e-05 ***
mean_elevation_scaled       -4.60340    1.81547     1.81571   2.535   0.0112 *  
number_of_checklists_scaled 20.82193    0.64287     0.64296  32.384  < 2e-16 ***
open_forest                 12.96359    0.98655     0.98669  13.138  < 2e-16 ***
permanent_water             -0.36087    0.59918     0.59923   0.602   0.5470    
realmAfrotropic              3.42908    1.91360     1.91387   1.792   0.0732 .  
realmAustralasia            -3.14597    3.13851     3.13896   1.002   0.3162    
realmIndomalayan            -1.55901    1.21340     1.21357   1.285   0.1989    
realmNearctic               -5.99756    1.41243     1.41263   4.246 2.18e-05 ***
realmNeotropic              -2.12324    1.42796     1.42817   1.487   0.1371    
shrubs                       6.29980    1.45822     1.45842   4.320 1.56e-05 ***
urban                        0.04047    0.18754     0.18756   0.216   0.8292    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  14.0932     0.8839      0.8840  15.942  < 2e-16 ***
closed_forest                 3.8989     0.7068      0.7069   5.516  < 2e-16 ***
cultivated                   11.3740     0.9813      0.9814  11.589  < 2e-16 ***
elevation_delta_scaled        0.9802     0.8683      0.8685   1.129   0.2591    
herbaceous_vegetation         9.8810     1.0889      1.0891   9.073  < 2e-16 ***
herbaceous_wetland           12.3517     3.1381      3.1386   3.935 8.30e-05 ***
mean_elevation_scaled        -4.6034     1.8155      1.8157   2.535   0.0112 *  
number_of_checklists_scaled  20.8219     0.6429      0.6430  32.384  < 2e-16 ***
open_forest                  12.9636     0.9866      0.9867  13.138  < 2e-16 ***
permanent_water              -0.6993     0.6775      0.6776   1.032   0.3021    
realmAfrotropic               3.4291     1.9136      1.9139   1.792   0.0732 .  
realmAustralasia             -3.1460     3.1385      3.1390   1.002   0.3162    
realmIndomalayan             -1.5590     1.2134      1.2136   1.285   0.1989    
realmNearctic                -5.9976     1.4124      1.4126   4.246 2.18e-05 ***
realmNeotropic               -2.1232     1.4280      1.4282   1.487   0.1371    
shrubs                        6.2998     1.4582      1.4584   4.320 1.56e-05 ***
urban                         0.1554     0.3423      0.3424   0.454   0.6499    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

merlin_mixed_model <- lmer(percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1|city_name), merlin)
r.squaredGLMM(merlin_mixed_model)

           R2m       R2c
[1,] 0.1614687 0.6330958

merlin_fixed_model <- lm(percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled, merlin)
anova(merlin_mixed_model, merlin_fixed_model)

refitting model(s) with ML (instead of REML)

Data: merlin
Models:
merlin_fixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled
merlin_mixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1 | city_name)
                   npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
merlin_fixed_model   17 50580 50700 -25273    50546                         
merlin_mixed_model   18 49305 49432 -24634    49269 1277.7  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(merlin_mixed_model, lm(percentage_of_regional_pool_present ~ city_name, merlin))

refitting model(s) with ML (instead of REML)

Data: merlin
Models:
merlin_mixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1 | city_name)
lm(percentage_of_regional_pool_present ~ city_name, merlin): percentage_of_regional_pool_present ~ city_name
                                                            npar   AIC   BIC logLik deviance Chisq  Df Pr(>Chisq)
merlin_mixed_model                                            18 49305 49432 -24634    49269                     
lm(percentage_of_regional_pool_present ~ city_name, merlin)  138 50605 51576 -25164    50329     0 120          1

merlin_richness_sum <- model_summary('merlin', 'richness', merlin_richness_result)
merlin_richness_sum

merlin_city_response <- create_city_dataset(merlin_mixed_model)

Joining, by = "name"

merlin_city_response[order(merlin_city_response$response), c('name', 'response')]

write_csv(merlin_city_response[, c('name', 'response')], 'percentage_of_regional_richness__output__merlin_city_richness_intercept.csv')

summary(merlin_mixed_model)

Linear mixed model fit by REML ['lmerMod']
Formula: percentage_of_regional_pool_present ~ mean_elevation_scaled +      elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation +  
    herbaceous_wetland + open_forest + permanent_water + shrubs +      realm + number_of_checklists_scaled + (1 | city_name)
   Data: merlin

REML criterion at convergence: 49236

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0484 -0.6453 -0.0591  0.5795  6.2777 

Random effects:
 Groups    Name        Variance Std.Dev.
 city_name (Intercept) 24.59    4.959   
 Residual              19.13    4.374   
Number of obs: 8443, groups:  city_name, 137

Fixed effects:
                            Estimate Std. Error t value
(Intercept)                  14.1438     0.8719  16.222
mean_elevation_scaled        -5.0263     1.8278  -2.750
elevation_delta_scaled        1.0054     0.8685   1.158
closed_forest                 3.8365     0.6942   5.526
cultivated                   11.3259     0.9672  11.710
herbaceous_vegetation         9.8109     1.0795   9.088
herbaceous_wetland           12.4718     3.1352   3.978
open_forest                  12.8555     0.9867  13.029
permanent_water              -0.7345     0.6589  -1.115
shrubs                        6.2175     1.4540   4.276
realmAfrotropic               3.5255     1.9121   1.844
realmAustralasia             -3.1666     3.1348  -1.010
realmIndomalayan             -1.5482     1.2120  -1.277
realmNearctic                -5.9691     1.4107  -4.231
realmNeotropic               -2.0731     1.4264  -1.453
number_of_checklists_scaled  20.8490     0.6433  32.411

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

merlin_richness_plot <- merlin_richness_sum[merlin_richness_sum$model == 'full',]
merlin_richness_plot$explanatory <- factor(merlin_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
merlin_richness_plot$type <- 'Proportion of landcover\n5 km around site'
merlin_richness_plot$type[merlin_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
merlin_richness_plot$type[merlin_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
merlin_richness_plot$type[merlin_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

ggplot(merlin_richness_plot, aes(y=explanatory, x=richness_merlin_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_merlin_estimate-richness_merlin_error, xmax=richness_merlin_estimate+richness_merlin_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(merlin_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")

ggsave('percentage_of_regional_pool__output__merlin_richness_result.jpg')

Saving 7.29 x 4.51 in image

Birdlife: Difference in richness

bl_richness_result <- model_average_richness(birdlife)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

bl_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <5 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/4/5/6/7/8/9/10/11      17 -25622.28 51278.63  0.00   0.32
1/2/3/4/5/6/7/8/9/10/11    18 -25621.28 51278.64  0.00   0.32
1/2/4/5/6/7/8/9/10/11/12   18 -25621.95 51279.97  1.34   0.16
1/2/3/4/5/6/7/8/9/10/11/12 19 -25620.96 51280.01  1.38   0.16
1/2/4/5/6/7/8/10/11/12     17 -25624.23 51282.53  3.89   0.05

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.3727     0.6656      0.6657  18.586  < 2e-16 ***
closed_forest                 4.6122     0.8136      0.8137   5.668  < 2e-16 ***
cultivated                   11.3601     1.1041      1.1043  10.287  < 2e-16 ***
herbaceous_vegetation        11.6039     1.2364      1.2366   9.384  < 2e-16 ***
herbaceous_wetland           12.0168     3.5407      3.5412   3.393  0.00069 ***
mean_elevation_scaled        -3.7033     1.5005      1.5007   2.468  0.01360 *  
number_of_checklists_scaled  25.0703     0.7297      0.7298  34.353  < 2e-16 ***
open_forest                  16.1070     1.1164      1.1165  14.426  < 2e-16 ***
permanent_water              -1.5753     0.8385      0.8386   1.879  0.06031 .  
realmAfrotropic              -1.4162     1.4153      1.4155   1.000  0.31709    
realmAustralasia             -3.2249     2.1301      2.1304   1.514  0.13009    
realmIndomalayan             -4.3603     0.8615      0.8616   5.061 4.00e-07 ***
realmNearctic                -2.6755     0.9471      0.9473   2.824  0.00474 ** 
realmNeotropic               -4.1732     1.0243      1.0245   4.074 4.63e-05 ***
shrubs                        9.0283     1.6084      1.6086   5.613  < 2e-16 ***
elevation_delta_scaled        0.2427     0.6903      0.6904   0.351  0.72523    
urban                         0.1323     0.3002      0.3002   0.441  0.65943    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.3727     0.6656      0.6657  18.586  < 2e-16 ***
closed_forest                 4.6122     0.8136      0.8137   5.668  < 2e-16 ***
cultivated                   11.3601     1.1041      1.1043  10.287  < 2e-16 ***
herbaceous_vegetation        11.6039     1.2364      1.2366   9.384  < 2e-16 ***
herbaceous_wetland           12.0168     3.5407      3.5412   3.393  0.00069 ***
mean_elevation_scaled        -3.7033     1.5005      1.5007   2.468  0.01360 *  
number_of_checklists_scaled  25.0703     0.7297      0.7298  34.353  < 2e-16 ***
open_forest                  16.1070     1.1164      1.1165  14.426  < 2e-16 ***
permanent_water              -1.6499     0.7830      0.7832   2.107  0.03514 *  
realmAfrotropic              -1.4162     1.4153      1.4155   1.000  0.31709    
realmAustralasia             -3.2249     2.1301      2.1304   1.514  0.13009    
realmIndomalayan             -4.3603     0.8615      0.8616   5.061 4.00e-07 ***
realmNearctic                -2.6755     0.9471      0.9473   2.824  0.00474 ** 
realmNeotropic               -4.1732     1.0243      1.0245   4.074 4.63e-05 ***
shrubs                        9.0283     1.6084      1.6086   5.613  < 2e-16 ***
elevation_delta_scaled        0.5102     0.9303      0.9304   0.548  0.58343    
urban                         0.3606     0.4040      0.4041   0.892  0.37217    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

birdlife_mixed_model <- lmer(percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + shrubs + realm + (1|city_name), birdlife)
r.squaredGLMM(birdlife_mixed_model)

           R2m       R2c
[1,] 0.1072897 0.3305271

summary(birdlife_mixed_model)

Linear mixed model fit by REML ['lmerMod']
Formula: percentage_of_regional_pool_present ~ mean_elevation_scaled +      closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland +  
    open_forest + shrubs + realm + (1 | city_name)
   Data: birdlife

REML criterion at convergence: 52351.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2243 -0.6782 -0.0908  0.6152  4.9342 

Random effects:
 Groups    Name        Variance Std.Dev.
 city_name (Intercept)  9.381   3.063   
 Residual              28.132   5.304   
Number of obs: 8443, groups:  city_name, 137

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            13.2793     0.6152  21.584
mean_elevation_scaled  -3.5910     1.4355  -2.502
closed_forest           3.6458     0.8310   4.387
cultivated             11.7687     1.1288  10.426
herbaceous_vegetation  10.1846     1.2886   7.904
herbaceous_wetland      9.8345     3.7552   2.619
open_forest            15.7218     1.1651  13.495
shrubs                  9.3424     1.6824   5.553
realmAfrotropic        -1.5413     1.4110  -1.092
realmAustralasia       -3.0458     2.0803  -1.464
realmIndomalayan       -4.2258     0.8522  -4.958
realmNearctic          -1.7182     0.9236  -1.860
realmNeotropic         -4.0139     1.0135  -3.960

Correlation matrix not shown by default, as p = 13 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

birdlife_fixed_model <- lm(percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + shrubs + realm, birdlife)
anova(birdlife_mixed_model, birdlife_fixed_model)

refitting model(s) with ML (instead of REML)

Data: birdlife
Models:
birdlife_fixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + shrubs + realm
birdlife_mixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + shrubs + realm + (1 | city_name)
                     npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
birdlife_fixed_model   14 53110 53208 -26541    53082                         
birdlife_mixed_model   15 52408 52514 -26189    52378 703.39  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(birdlife_fixed_model)

Call:
lm(formula = percentage_of_regional_pool_present ~ mean_elevation_scaled + 
    closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + 
    open_forest + shrubs + realm, data = birdlife)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.7661  -4.0084  -0.5982   3.3769  27.5192 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            14.5614     0.2319  62.795  < 2e-16 ***
mean_elevation_scaled  -3.7600     0.5323  -7.064 1.75e-12 ***
closed_forest           2.0352     0.7227   2.816  0.00487 ** 
cultivated             16.0866     0.9808  16.401  < 2e-16 ***
herbaceous_vegetation   9.2569     1.0864   8.521  < 2e-16 ***
herbaceous_wetland      0.2377     3.5321   0.067  0.94636    
open_forest            14.3239     1.0159  14.100  < 2e-16 ***
shrubs                  6.7982     1.1102   6.124 9.56e-10 ***
realmAfrotropic        -2.1892     0.7785  -2.812  0.00493 ** 
realmAustralasia       -4.3652     0.3482 -12.537  < 2e-16 ***
realmIndomalayan       -5.7999     0.3142 -18.460  < 2e-16 ***
realmNearctic          -2.5121     0.2285 -10.995  < 2e-16 ***
realmNeotropic         -5.2385     0.3279 -15.978  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.615 on 8430 degrees of freedom
Multiple R-squared:  0.1424,    Adjusted R-squared:  0.1411 
F-statistic: 116.6 on 12 and 8430 DF,  p-value: < 2.2e-16

birdlife_city_response <- create_city_dataset(birdlife_mixed_model)

Joining, by = "name"

birdlife_city_response[order(birdlife_city_response$response), c('name', 'response')]

write_csv(birdlife_city_response[, c('name', 'response')], 'percentage_of_regional_richness__output__birdlife_city_richness_intercept.csv')

bl_richness_sum <- model_summary('birdlife', 'richness', bl_richness_result)
bl_richness_sum

summary(birdlife_mixed_model)

Linear mixed model fit by REML ['lmerMod']
Formula: percentage_of_regional_pool_present ~ mean_elevation_scaled +      closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland +  
    open_forest + shrubs + realm + (1 | city_name)
   Data: birdlife

REML criterion at convergence: 52351.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2243 -0.6782 -0.0908  0.6152  4.9342 

Random effects:
 Groups    Name        Variance Std.Dev.
 city_name (Intercept)  9.381   3.063   
 Residual              28.132   5.304   
Number of obs: 8443, groups:  city_name, 137

Fixed effects:
                      Estimate Std. Error t value
(Intercept)            13.2793     0.6152  21.584
mean_elevation_scaled  -3.5910     1.4355  -2.502
closed_forest           3.6458     0.8310   4.387
cultivated             11.7687     1.1288  10.426
herbaceous_vegetation  10.1846     1.2886   7.904
herbaceous_wetland      9.8345     3.7552   2.619
open_forest            15.7218     1.1651  13.495
shrubs                  9.3424     1.6824   5.553
realmAfrotropic        -1.5413     1.4110  -1.092
realmAustralasia       -3.0458     2.0803  -1.464
realmIndomalayan       -4.2258     0.8522  -4.958
realmNearctic          -1.7182     0.9236  -1.860
realmNeotropic         -4.0139     1.0135  -3.960

Correlation matrix not shown by default, as p = 13 > 12.
Use print(x, correlation=TRUE)  or
    vcov(x)        if you need it

birdlife_richness_plot <- bl_richness_sum[bl_richness_sum$model == 'full',]
birdlife_richness_plot$explanatory <- factor(birdlife_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
birdlife_richness_plot$type <- 'Proportion of landcover\n5 km around site'
birdlife_richness_plot$type[birdlife_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
birdlife_richness_plot$type[birdlife_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
birdlife_richness_plot$type[birdlife_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

birdlife_plot = ggplot(birdlife_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(birdlife_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")
birdlife_plot
ggsave('percentage_of_regional_pool__output__birdlife_richness_result.jpg')

Saving 7.29 x 4.51 in image

Both: Difference in richness

both_richness_result <- model_average_richness(both)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"
Warning: Model failed to converge with max|grad| = 0.00314894 (tol = 0.002, component 1)

both_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <7 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/3/4/5/6/7/8/9/10/11    18 -26535.43 53106.94  0.00   0.35
1/2/4/5/6/7/8/9/10/11      17 -26536.96 53107.99  1.06   0.20
1/2/3/4/5/6/7/8/9/10/11/12 19 -26535.19 53108.46  1.53   0.16
1/2/4/5/6/7/8/9/10/11/12   18 -26536.70 53109.48  2.54   0.10
1/2/3/4/5/6/7/8/10/11      17 -26537.95 53109.98  3.05   0.08
1/2/3/4/5/6/7/8/10/11/12   18 -26537.13 53110.34  3.41   0.06
1/2/4/5/6/7/8/10/11        16 -26539.40 53110.85  3.92   0.05

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 16.76949    0.93298     0.93311  17.972  < 2e-16 ***
closed_forest                5.23725    0.89790     0.89803   5.832  < 2e-16 ***
cultivated                  13.56656    1.23804     1.23822  10.957  < 2e-16 ***
elevation_delta_scaled       0.70796    1.00687     1.00697   0.703  0.48202    
herbaceous_vegetation       12.60120    1.37494     1.37514   9.164  < 2e-16 ***
herbaceous_wetland          16.36089    3.95058     3.95114   4.141 3.46e-05 ***
mean_elevation_scaled       -5.47279    2.05661     2.05689   2.661  0.00780 ** 
number_of_checklists_scaled 26.72527    0.81075     0.81087  32.959  < 2e-16 ***
open_forest                 17.44821    1.24768     1.24785  13.983  < 2e-16 ***
permanent_water             -1.20865    0.97092     0.97101   1.245  0.21323    
realmAfrotropic              1.82107    2.01065     2.01093   0.906  0.36516    
realmAustralasia            -0.58699    3.20839     3.20885   0.183  0.85485    
realmIndomalayan            -3.30686    1.25716     1.25734   2.630  0.00854 ** 
realmNearctic               -6.40401    1.43543     1.43563   4.461 8.20e-06 ***
realmNeotropic              -4.08496    1.48660     1.48681   2.747  0.00601 ** 
shrubs                       9.75278    1.82036     1.82062   5.357 1.00e-07 ***
urban                        0.09439    0.28888     0.28892   0.327  0.74389    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  16.7695     0.9330      0.9331  17.972  < 2e-16 ***
closed_forest                 5.2372     0.8979      0.8980   5.832  < 2e-16 ***
cultivated                   13.5666     1.2380      1.2382  10.957  < 2e-16 ***
elevation_delta_scaled        1.0912     1.0698      1.0699   1.020  0.30776    
herbaceous_vegetation        12.6012     1.3749      1.3751   9.164  < 2e-16 ***
herbaceous_wetland           16.3609     3.9506      3.9511   4.141 3.46e-05 ***
mean_elevation_scaled        -5.4728     2.0566      2.0569   2.661  0.00780 ** 
number_of_checklists_scaled  26.7253     0.8108      0.8109  32.959  < 2e-16 ***
open_forest                  17.4482     1.2477      1.2479  13.983  < 2e-16 ***
permanent_water              -1.4886     0.8628      0.8629   1.725  0.08451 .  
realmAfrotropic               1.8211     2.0106      2.0109   0.906  0.36516    
realmAustralasia             -0.5870     3.2084      3.2089   0.183  0.85485    
realmIndomalayan             -3.3069     1.2572      1.2573   2.630  0.00854 ** 
realmNearctic                -6.4040     1.4354      1.4356   4.461 8.20e-06 ***
realmNeotropic               -4.0850     1.4866      1.4868   2.747  0.00601 ** 
shrubs                        9.7528     1.8204      1.8206   5.357 1.00e-07 ***
urban                         0.2924     0.4479      0.4480   0.653  0.51389    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

both_mixed_model <- lmer(percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1|city_name), both)
r.squaredGLMM(both_mixed_model)

           R2m       R2c
[1,] 0.1688506 0.5417001

both_fixed_model <- lm(percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled, both)
anova(both_mixed_model, both_fixed_model)

refitting model(s) with ML (instead of REML)

Data: both
Models:
both_fixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled
both_mixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + elevation_delta_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1 | city_name)
                 npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
both_fixed_model   17 54137 54257 -27052    54103                         
both_mixed_model   18 53144 53271 -26554    53108 994.91  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

both_richness_sum <- model_summary('both', 'richness', both_richness_result)
both_richness_sum

both_city_response <- create_city_dataset(both_mixed_model)

Joining, by = "name"

both_city_response[order(both_city_response$response), c('name', 'response')]

write_csv(both_city_response[, c('name', 'response')], 'percentage_of_regional_richness__output__both_city_richness_intercept.csv')

Either: Difference in richness

either_richness_result <- model_average_richness(either)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

either_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <8 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/4/5/6/7/8/9/10/11      17 -23755.12 47544.31  0.00   0.23
1/2/3/4/5/6/7/8/9/10/11    18 -23754.19 47544.46  0.15   0.22
1/2/4/5/6/7/8/10/11        16 -23756.62 47545.30  0.98   0.14
1/2/3/4/5/6/7/8/10/11      17 -23755.74 47545.55  1.24   0.13
1/2/4/5/6/7/8/9/10/11/12   18 -23755.19 47546.46  2.15   0.08
1/2/3/4/5/6/7/8/9/10/11/12 19 -23754.28 47546.64  2.33   0.07
1/2/4/5/6/7/8/10/11/12     17 -23756.30 47546.67  2.35   0.07
1/2/3/4/5/6/7/8/10/11/12   18 -23755.43 47546.95  2.63   0.06

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 10.79309    0.56100     0.56108  19.236  < 2e-16 ***
closed_forest                3.48703    0.64338     0.64347   5.419  1.0e-07 ***
cultivated                   9.95937    0.87880     0.87892  11.331  < 2e-16 ***
herbaceous_vegetation        8.97326    0.98504     0.98518   9.108  < 2e-16 ***
herbaceous_wetland          10.00181    2.83810     2.83851   3.524 0.000426 ***
mean_elevation_scaled       -3.43498    1.27445     1.27463   2.695 0.007041 ** 
number_of_checklists_scaled 19.74517    0.58407     0.58416  33.801  < 2e-16 ***
open_forest                 12.13158    0.89205     0.89218  13.598  < 2e-16 ***
permanent_water             -0.51126    0.63442     0.63447   0.806 0.420351    
realmAfrotropic              0.07749    1.21551     1.21568   0.064 0.949174    
realmAustralasia            -1.64760    1.86786     1.86813   0.882 0.377802    
realmIndomalayan            -2.89610    0.74686     0.74697   3.877 0.000106 ***
realmNearctic               -3.08571    0.83146     0.83158   3.711 0.000207 ***
realmNeotropic              -2.67969    0.88656     0.88669   3.022 0.002510 ** 
shrubs                       6.10974    1.29325     1.29344   4.724  2.3e-06 ***
elevation_delta_scaled       0.26239    0.58738     0.58745   0.447 0.655126    
urban                        0.07007    0.20080     0.20082   0.349 0.727156    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 10.79309    0.56100     0.56108  19.236  < 2e-16 ***
closed_forest                3.48703    0.64338     0.64347   5.419  1.0e-07 ***
cultivated                   9.95937    0.87880     0.87892  11.331  < 2e-16 ***
herbaceous_vegetation        8.97326    0.98504     0.98518   9.108  < 2e-16 ***
herbaceous_wetland          10.00181    2.83810     2.83851   3.524 0.000426 ***
mean_elevation_scaled       -3.43498    1.27445     1.27463   2.695 0.007041 ** 
number_of_checklists_scaled 19.74517    0.58407     0.58416  33.801  < 2e-16 ***
open_forest                 12.13158    0.89205     0.89218  13.598  < 2e-16 ***
permanent_water             -0.85325    0.61637     0.61645   1.384 0.166319    
realmAfrotropic              0.07749    1.21551     1.21568   0.064 0.949174    
realmAustralasia            -1.64760    1.86786     1.86813   0.882 0.377802    
realmIndomalayan            -2.89610    0.74686     0.74697   3.877 0.000106 ***
realmNearctic               -3.08571    0.83146     0.83158   3.711 0.000207 ***
realmNeotropic              -2.67969    0.88656     0.88669   3.022 0.002510 ** 
shrubs                       6.10974    1.29325     1.29344   4.724  2.3e-06 ***
elevation_delta_scaled       0.55247    0.75246     0.75257   0.734 0.462881    
urban                        0.24538    0.31334     0.31339   0.783 0.433635    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

either_mixed_model <- lmer(percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1|city_name), either)
r.squaredGLMM(either_mixed_model)

           R2m       R2c
[1,] 0.1617725 0.4424337

either_fixed_model <- lm(percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled, either)
anova(either_mixed_model, either_fixed_model)

refitting model(s) with ML (instead of REML)

Data: either
Models:
either_fixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled
either_mixed_model: percentage_of_regional_pool_present ~ mean_elevation_scaled + closed_forest + cultivated + herbaceous_vegetation + herbaceous_wetland + open_forest + permanent_water + shrubs + realm + number_of_checklists_scaled + (1 | city_name)
                   npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
either_fixed_model   16 48390 48503 -24179    48358                         
either_mixed_model   17 47567 47687 -23767    47533 824.93  1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

either_richness_sum <- model_summary('either', 'richness', either_richness_result)
either_richness_sum

either_city_response <- create_city_dataset(either_mixed_model)

Joining, by = "name"

either_city_response[order(either_city_response$response), c('name', 'response')]

write_csv(either_city_response[, c('name', 'response')], 'percentage_of_regional_richness__output__either_city_richness_intercept.csv')

Full result

richness_all <- full_join(full_join(merlin_richness_sum, bl_richness_sum), full_join(both_richness_sum, either_richness_sum))

Joining, by = c("explanatory", "model")Joining, by = c("explanatory", "model")Joining, by = c("explanatory", "model")

write_csv(richness_all, "percentage_of_regional_richness__output__richness_result.csv")
richness_all

Additional figures

library(tidyr)
library(ggplot2)
library(ggpubr)

environment <- merlin[,c('closed_forest', 'cultivated', 'herbaceous_vegetation', 'herbaceous_wetland', 'open_forest', 'permanent_water', 'shrubs', 'urban')]
names(environment) <- c('Closed Forest', 'Cultivated', 'Herbaceous Vegetation', 'Herbaceous Wetland', 'Open Forest', 'Permanent Water', 'Shurbs', 'Urban')
landcover <- gather(environment, landcover, proportion)

ggplot(landcover, aes(x = landcover, y = proportion)) + geom_boxplot() + xlab('Landcover Class') + ylab('Proportion in 5 km around sites') +
  scale_x_discrete(labels = scales::wrap_format(10)) + theme_bw() + labs(title = "Landcover Proportions")

ggsave("percentage_of_regional_richness__output__som_landcover_proportions.jpg")

Saving 7.29 x 4.51 in image

ggarrange(
  ggplot(merlin, aes(y = elevation_delta)) + geom_boxplot() + xlab('Elevation Delta') + ylab('Average Metres in 5 km around sites') +
    scale_x_discrete(labels = scales::wrap_format(10)) + theme_bw(),
  ggplot(merlin, aes(y = mean_elevation)) + geom_boxplot() + xlab('Mean Elevation') + ylab('Average Metres in 5 km around sites') +
    scale_x_discrete(labels = scales::wrap_format(10)) + theme_bw()
)

ggsave("percentage_of_regional_richness__output__som_elevations.jpg")

Saving 7.29 x 4.51 in image

Birdlife [300 or less checklists]: Difference in richness

birdlife300 = birdlife[birdlife$number_of_checklists <= 300,]
nrow(birdlife300)

[1] 7190

blr300_richness_result <- model_average_richness(birdlife300)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

blr300_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <6 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/3/4/5/6/7/8/9/10/11/12 19 -21709.04 43456.18  0.00   0.25
1/2/4/5/6/7/8/9/10/11/12   18 -21710.05 43456.19  0.02   0.25
1/2/3/4/5/6/7/8/9/10/11    18 -21710.39 43456.88  0.70   0.18
1/2/4/5/6/7/8/9/10/11      17 -21711.44 43456.97  0.80   0.17
1/2/4/5/6/7/8/10/11/12     17 -21712.27 43458.62  2.44   0.07
1/2/3/4/5/6/7/8/10/11/12   18 -21711.30 43458.69  2.52   0.07

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  11.7790     0.7302      0.7304  16.128  < 2e-16 ***
closed_forest                 4.7008     0.8989      0.8990   5.229  2.0e-07 ***
cultivated                   12.3055     1.1979      1.1980  10.271  < 2e-16 ***
elevation_delta_scaled        0.2506     0.7329      0.7330   0.342 0.732465    
herbaceous_vegetation        11.5511     1.3443      1.3445   8.591  < 2e-16 ***
herbaceous_wetland           13.4647     3.7010      3.7016   3.638 0.000275 ***
mean_elevation_scaled        -3.4077     1.5165      1.5168   2.247 0.024661 *  
number_of_checklists_scaled  29.3514     1.9470      1.9473  15.073  < 2e-16 ***
open_forest                  18.2162     1.1871      1.1873  15.342  < 2e-16 ***
permanent_water              -1.5377     1.0807      1.0808   1.423 0.154819    
realmAfrotropic              -1.3731     1.3897      1.3900   0.988 0.323228    
realmAustralasia             -3.1636     2.0863      2.0866   1.516 0.129492    
realmIndomalayan             -4.2550     0.8470      0.8471   5.023  5.0e-07 ***
realmNearctic                -2.7905     0.9289      0.9291   3.003 0.002669 ** 
realmNeotropic               -4.1342     1.0078      1.0079   4.102  4.1e-05 ***
shrubs                        8.6184     1.6968      1.6971   5.078  4.0e-07 ***
urban                         0.4908     0.5013      0.5013   0.979 0.327582    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  11.7790     0.7302      0.7304  16.128  < 2e-16 ***
closed_forest                 4.7008     0.8989      0.8990   5.229  2.0e-07 ***
cultivated                   12.3055     1.1979      1.1980  10.271  < 2e-16 ***
elevation_delta_scaled        0.4973     0.9713      0.9714   0.512 0.608703    
herbaceous_vegetation        11.5511     1.3443      1.3445   8.591  < 2e-16 ***
herbaceous_wetland           13.4647     3.7010      3.7016   3.638 0.000275 ***
mean_elevation_scaled        -3.4077     1.5165      1.5168   2.247 0.024661 *  
number_of_checklists_scaled  29.3514     1.9470      1.9473  15.073  < 2e-16 ***
open_forest                  18.2162     1.1871      1.1873  15.342  < 2e-16 ***
permanent_water              -1.8020     0.9447      0.9448   1.907 0.056488 .  
realmAfrotropic              -1.3731     1.3897      1.3900   0.988 0.323228    
realmAustralasia             -3.1636     2.0863      2.0866   1.516 0.129492    
realmIndomalayan             -4.2550     0.8470      0.8471   5.023  5.0e-07 ***
realmNearctic                -2.7905     0.9289      0.9291   3.003 0.002669 ** 
realmNeotropic               -4.1342     1.0078      1.0079   4.102  4.1e-05 ***
shrubs                        8.6184     1.6968      1.6971   5.078  4.0e-07 ***
urban                         0.7533     0.4335      0.4335   1.738 0.082294 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

bl300_richness_sum <- model_summary('birdlife', 'richness', blr300_richness_result)
bl300_richness_sum

birdlife300_richness_plot <- bl300_richness_sum[bl300_richness_sum$model == 'full',]

birdlife300_richness_plot$explanatory <- factor(birdlife300_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
birdlife300_richness_plot$type <- 'Proportion of landcover\n5 km around site'
birdlife300_richness_plot$type[birdlife300_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
birdlife300_richness_plot$type[birdlife300_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
birdlife300_richness_plot$type[birdlife300_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

birdlife300_plot = ggplot(birdlife300_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(birdlife300_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")
birdlife300_plot
ggsave('percentage_of_regional_richness__output__birdlife_richness_300_checklists_or_less_result.jpg')

Saving 7.29 x 4.51 in image

Birdlife [No Australasia]: Difference in richness

birdlifeNoAustralia = birdlife[birdlife$realm != 'Australasia',]
nrow(birdlifeNoAustralia)

[1] 8014

blrNoAustralia_richness_result <- model_average_richness(birdlifeNoAustralia)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

blrNoAustralia_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <4 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/4/5/6/7/8/9/10/11      16 -24381.60 48795.26  0.00   0.36
1/2/3/4/5/6/7/8/9/10/11    17 -24380.71 48795.49  0.23   0.32
1/2/4/5/6/7/8/9/10/11/12   17 -24381.36 48796.80  1.54   0.17
1/2/3/4/5/6/7/8/9/10/11/12 18 -24380.48 48797.04  1.78   0.15

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 12.47596    0.65663     0.65672  18.997  < 2e-16 ***
closed_forest                4.63809    0.82943     0.82955   5.591  < 2e-16 ***
cultivated                  11.60344    1.10345     1.10361  10.514  < 2e-16 ***
herbaceous_vegetation       10.41195    1.27598     1.27618   8.159  < 2e-16 ***
herbaceous_wetland          11.44994    3.58957     3.59011   3.189  0.00143 ** 
mean_elevation_scaled       -3.72351    1.50342     1.50364   2.476  0.01327 *  
number_of_checklists_scaled 24.96892    0.74474     0.74485  33.522  < 2e-16 ***
open_forest                 15.51587    1.13841     1.13858  13.627  < 2e-16 ***
permanent_water             -1.62990    0.78561     0.78572   2.074  0.03804 *  
realmAfrotropic             -1.37071    1.42023     1.42044   0.965  0.33455    
realmIndomalayan            -4.40896    0.86422     0.86435   5.101 3.00e-07 ***
realmNearctic               -2.63636    0.94943     0.94957   2.776  0.00550 ** 
realmNeotropic              -4.11584    1.02769     1.02785   4.004 6.22e-05 ***
shrubs                       9.54960    1.61915     1.61939   5.897  < 2e-16 ***
elevation_delta_scaled       0.10597    0.65561     0.65571   0.162  0.87161    
urban                        0.08344    0.26011     0.26014   0.321  0.74841    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.4760     0.6566      0.6567  18.997  < 2e-16 ***
closed_forest                 4.6381     0.8294      0.8296   5.591  < 2e-16 ***
cultivated                   11.6034     1.1034      1.1036  10.514  < 2e-16 ***
herbaceous_vegetation        10.4120     1.2760      1.2762   8.159  < 2e-16 ***
herbaceous_wetland           11.4499     3.5896      3.5901   3.189  0.00143 ** 
mean_elevation_scaled        -3.7235     1.5034      1.5036   2.476  0.01327 *  
number_of_checklists_scaled  24.9689     0.7447      0.7449  33.522  < 2e-16 ***
open_forest                  15.5159     1.1384      1.1386  13.627  < 2e-16 ***
permanent_water              -1.6299     0.7856      0.7857   2.074  0.03804 *  
realmAfrotropic              -1.3707     1.4202      1.4204   0.965  0.33455    
realmIndomalayan             -4.4090     0.8642      0.8644   5.101 3.00e-07 ***
realmNearctic                -2.6364     0.9494      0.9496   2.776  0.00550 ** 
realmNeotropic               -4.1158     1.0277      1.0278   4.004 6.22e-05 ***
shrubs                        9.5496     1.6191      1.6194   5.897  < 2e-16 ***
elevation_delta_scaled        0.2251     0.9414      0.9416   0.239  0.81103    
urban                         0.2641     0.4080      0.4080   0.647  0.51748    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

blNoAustralia_richness_sum <- model_summary('birdlife', 'richness', blrNoAustralia_richness_result)
blNoAustralia_richness_sum

birdlifeNoAustralia_richness_plot <- blNoAustralia_richness_sum[blNoAustralia_richness_sum$model == 'full',]

birdlifeNoAustralia_richness_plot$explanatory <- factor(birdlifeNoAustralia_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
birdlifeNoAustralia_richness_plot$type <- 'Proportion of landcover\n5 km around site'
birdlifeNoAustralia_richness_plot$type[birdlifeNoAustralia_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
birdlifeNoAustralia_richness_plot$type[birdlifeNoAustralia_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
birdlifeNoAustralia_richness_plot$type[birdlifeNoAustralia_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

birdlifeNoAustralia_plot = ggplot(birdlifeNoAustralia_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(birdlifeNoAustralia_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")
birdlifeNoAustralia_plot
ggsave('percentage_of_regional_richness__output__birdlife_richness_no_australia.jpg')

Saving 7.29 x 4.51 in image

Birdlife [No Afrotropics]: Difference in richness

birdlifeNoAfrica = birdlife[birdlife$realm != 'Afrotropic',]
nrow(birdlifeNoAfrica)

[1] 8385

blrNoAfrica_richness_result <- model_average_richness(birdlifeNoAfrica)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

blrNoAfrica_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <4 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/4/5/6/7/8/9/10/11      16 -25422.24 50876.54  0.00   0.36
1/2/3/4/5/6/7/8/9/10/11    17 -25421.29 50876.65  0.10   0.34
1/2/4/5/6/7/8/9/10/11/12   17 -25422.07 50878.22  1.67   0.16
1/2/3/4/5/6/7/8/9/10/11/12 18 -25421.13 50878.34  1.79   0.15

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 12.39958    0.66378     0.66387  18.678  < 2e-16 ***
closed_forest                4.56883    0.80025     0.80037   5.708  < 2e-16 ***
cultivated                  11.11759    1.09826     1.09841  10.121  < 2e-16 ***
herbaceous_vegetation       11.98855    1.24148     1.24165   9.655  < 2e-16 ***
herbaceous_wetland          11.07292    3.54284     3.54335   3.125  0.00178 ** 
mean_elevation_scaled       -3.27277    1.60646     1.60669   2.037  0.04165 *  
number_of_checklists_scaled 25.05393    0.72773     0.72784  34.422  < 2e-16 ***
open_forest                 16.04934    1.11065     1.11081  14.448  < 2e-16 ***
permanent_water             -1.64144    0.77354     0.77365   2.122  0.03386 *  
realmAustralasia            -3.17949    2.18437     2.18469   1.455  0.14557    
realmIndomalayan            -4.32224    0.88038     0.88051   4.909 9.00e-07 ***
realmNearctic               -2.68378    0.97151     0.97165   2.762  0.00574 ** 
realmNeotropic              -4.26194    1.05094     1.05110   4.055 5.02e-05 ***
shrubs                       8.49079    1.61473     1.61496   5.258 1.00e-07 ***
elevation_delta_scaled       0.19823    0.68537     0.68546   0.289  0.77243    
urban                        0.06957    0.24221     0.24224   0.287  0.77398    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.3996     0.6638      0.6639  18.678  < 2e-16 ***
closed_forest                 4.5688     0.8003      0.8004   5.708  < 2e-16 ***
cultivated                   11.1176     1.0983      1.0984  10.121  < 2e-16 ***
herbaceous_vegetation        11.9886     1.2415      1.2417   9.655  < 2e-16 ***
herbaceous_wetland           11.0729     3.5428      3.5434   3.125  0.00178 ** 
mean_elevation_scaled        -3.2728     1.6065      1.6067   2.037  0.04165 *  
number_of_checklists_scaled  25.0539     0.7277      0.7278  34.422  < 2e-16 ***
open_forest                  16.0493     1.1106      1.1108  14.448  < 2e-16 ***
permanent_water              -1.6414     0.7735      0.7737   2.122  0.03386 *  
realmAustralasia             -3.1795     2.1844      2.1847   1.455  0.14557    
realmIndomalayan             -4.3222     0.8804      0.8805   4.909 9.00e-07 ***
realmNearctic                -2.6838     0.9715      0.9717   2.762  0.00574 ** 
realmNeotropic               -4.2619     1.0509      1.0511   4.055 5.02e-05 ***
shrubs                        8.4908     1.6147      1.6150   5.258 1.00e-07 ***
elevation_delta_scaled        0.4076     0.9383      0.9385   0.434  0.66406    
urban                         0.2308     0.3968      0.3968   0.582  0.56085    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

blNoAfrica_richness_sum = model_summary('birdlife', 'richness', blrNoAfrica_richness_result)

birdlifeNoAfrica_richness_plot <- blNoAfrica_richness_sum[blNoAfrica_richness_sum$model == 'full',]

birdlifeNoAfrica_richness_plot$explanatory <- factor(birdlifeNoAfrica_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
birdlifeNoAfrica_richness_plot$type <- 'Proportion of landcover\n5 km around site'
birdlifeNoAfrica_richness_plot$type[birdlifeNoAfrica_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
birdlifeNoAfrica_richness_plot$type[birdlifeNoAfrica_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
birdlifeNoAfrica_richness_plot$type[birdlifeNoAfrica_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

birdlifeNoAfrica_plot = ggplot(birdlifeNoAfrica_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(birdlifeNoAfrica_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")
birdlifeNoAfrica_plot
ggsave('percentage_of_regional_richness__output__birdlife_richness_no_africa.jpg')

Saving 7.29 x 4.51 in image

Birdlife [No Afrotropics and no Australisia]: Difference in richness

birdlifeNeitherAfrNorAus = birdlife[birdlife$realm != 'Afrotropic' & birdlife$realm != 'Australasia',]
nrow(birdlifeNeitherAfrNorAus)

[1] 7956

blrNAA_richness_result <- model_average_richness(birdlifeNeitherAfrNorAus)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

blrNAA_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <4 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/4/5/6/7/8/9/10/11      15 -24182.21 48394.49  0.00   0.38
1/2/3/4/5/6/7/8/9/10/11    16 -24181.34 48394.74  0.26   0.33
1/2/4/5/6/7/8/9/10/11/12   16 -24182.11 48396.28  1.79   0.15
1/2/3/4/5/6/7/8/9/10/11/12 17 -24181.23 48396.54  2.05   0.14

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                 12.47873    0.66382     0.66392  18.795  < 2e-16 ***
closed_forest                4.61632    0.82369     0.82382   5.604  < 2e-16 ***
cultivated                  11.39742    1.10640     1.10657  10.300  < 2e-16 ***
herbaceous_vegetation       10.81810    1.28914     1.28933   8.390  < 2e-16 ***
herbaceous_wetland          10.50456    3.59220     3.59275   2.924  0.00346 ** 
mean_elevation_scaled       -3.30266    1.61097     1.61122   2.050  0.04038 *  
number_of_checklists_scaled 24.94747    0.74290     0.74301  33.576  < 2e-16 ***
open_forest                 15.48966    1.13689     1.13707  13.622  < 2e-16 ***
permanent_water             -1.61680    0.77844     0.77855   2.077  0.03783 *  
realmIndomalayan            -4.37384    0.88300     0.88313   4.953 7.00e-07 ***
realmNearctic               -2.64602    0.97363     0.97378   2.717  0.00658 ** 
realmNeotropic              -4.20152    1.05432     1.05448   3.984 6.76e-05 ***
shrubs                       9.02543    1.62801     1.62826   5.543  < 2e-16 ***
elevation_delta_scaled       0.05651    0.65242     0.65252   0.087  0.93099    
urban                        0.04821    0.23246     0.23249   0.207  0.83571    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.4787     0.6638      0.6639  18.795  < 2e-16 ***
closed_forest                 4.6163     0.8237      0.8238   5.604  < 2e-16 ***
cultivated                   11.3974     1.1064      1.1066  10.300  < 2e-16 ***
herbaceous_vegetation        10.8181     1.2891      1.2893   8.390  < 2e-16 ***
herbaceous_wetland           10.5046     3.5922      3.5928   2.924  0.00346 ** 
mean_elevation_scaled        -3.3027     1.6110      1.6112   2.050  0.04038 *  
number_of_checklists_scaled  24.9475     0.7429      0.7430  33.576  < 2e-16 ***
open_forest                  15.4897     1.1369      1.1371  13.622  < 2e-16 ***
permanent_water              -1.6168     0.7784      0.7786   2.077  0.03783 *  
realmIndomalayan             -4.3738     0.8830      0.8831   4.953 7.00e-07 ***
realmNearctic                -2.6460     0.9736      0.9738   2.717  0.00658 ** 
realmNeotropic               -4.2015     1.0543      1.0545   3.984 6.76e-05 ***
shrubs                        9.0254     1.6280      1.6283   5.543  < 2e-16 ***
elevation_delta_scaled        0.1208     0.9498      0.9499   0.127  0.89882    
urban                         0.1665     0.4085      0.4086   0.407  0.68367    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

blNAA_richness_sum <- model_summary('birdlife', 'richness', blrNAA_richness_result)
blNAA_richness_sum

birdlifeNAA_richness_plot <- blNAA_richness_sum[blNAA_richness_sum$model == 'full',]

birdlifeNAA_richness_plot$explanatory <- factor(birdlifeNAA_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
birdlifeNAA_richness_plot$type <- 'Proportion of landcover\n5 km around site'
birdlifeNAA_richness_plot$type[birdlifeNAA_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
birdlifeNAA_richness_plot$type[birdlifeNAA_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
birdlifeNAA_richness_plot$type[birdlifeNAA_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

birdlifeNAA_plot = ggplot(birdlifeNAA_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(birdlifeNAA_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")
birdlifeNAA_plot
ggsave('percentage_of_regional_richness__output__birdlife_richness_no_africa_nor_australisia.jpg')

Saving 7.29 x 4.51 in image

Birdlife [Urban Only]: Difference in richness

birdlife$urban_rounded = floor(birdlife$urban * 10) / 10
ggplot(birdlife, aes(x = urban_rounded)) + geom_bar()

birdlifeHigherUrban = birdlife[birdlife$urban > 0.2,]
nrow(birdlifeHigherUrban)

[1] 8262

blrHigherUrban_richness_result <- model_average_richness(birdlifeHigherUrban)

Warning: comparing models fitted by REMLFixed term is "(Intercept)"
Warning: Model failed to converge with max|grad| = 0.00435943 (tol = 0.002, component 1)

blrHigherUrban_richness_result

Call:
model.avg(object = dredge_result, subset = delta < aic_delta)

Component model call: 
lmer(formula = percentage_of_regional_pool_present ~ <8 unique rhs>, data = data)

Component models: 
                           df    logLik     AICc delta weight
1/2/3/4/5/6/7/8/9/10/11    18 -25059.71 50155.49  0.00   0.19
1/2/4/5/6/7/8/9/10/11      17 -25060.76 50155.59  0.10   0.18
1/2/3/4/5/6/7/8/9/10/11/12 19 -25059.12 50156.34  0.84   0.12
1/2/4/5/6/7/8/9/10/11/12   18 -25060.15 50156.39  0.90   0.12
1/2/4/5/6/7/8/10/11/12     17 -25061.29 50156.65  1.16   0.11
1/2/3/4/5/6/7/8/10/11/12   18 -25060.29 50156.65  1.16   0.11
1/2/3/4/5/6/7/8/10/11      17 -25061.47 50157.02  1.53   0.09
1/2/4/5/6/7/8/10/11        16 -25062.49 50157.05  1.55   0.09

Term codes: 
              closed_forest                  cultivated      elevation_delta_scaled       herbaceous_vegetation          herbaceous_wetland       mean_elevation_scaled 
                          1                           2                           3                           4                           5                           6 
number_of_checklists_scaled                 open_forest             permanent_water                       realm                      shrubs                       urban 
                          7                           8                           9                          10                          11                          12 

Model-averaged coefficients:  
(full average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.0109     0.6931      0.6931  17.328  < 2e-16 ***
closed_forest                 4.8906     0.8534      0.8535   5.730  < 2e-16 ***
cultivated                   12.1076     1.1800      1.1801  10.259  < 2e-16 ***
elevation_delta_scaled        0.2824     0.7237      0.7238   0.390 0.696461    
herbaceous_vegetation        11.4514     1.2644      1.2646   9.055  < 2e-16 ***
herbaceous_wetland           11.7865     3.6101      3.6106   3.264 0.001097 ** 
mean_elevation_scaled        -3.9688     1.5221      1.5223   2.607 0.009134 ** 
number_of_checklists_scaled  24.8715     0.7373      0.7374  33.729  < 2e-16 ***
open_forest                  17.4451     1.1501      1.1502  15.167  < 2e-16 ***
permanent_water              -0.6019     0.8028      0.8029   0.750 0.453499    
realmAfrotropic              -0.6229     1.4407      1.4409   0.432 0.665529    
realmAustralasia             -3.0481     2.1240      2.1243   1.435 0.151329    
realmIndomalayan             -4.2920     0.8641      0.8642   4.967    7e-07 ***
realmNearctic                -2.5242     0.9463      0.9464   2.667 0.007648 ** 
realmNeotropic               -3.9439     1.0244      1.0245   3.849 0.000118 ***
shrubs                        9.1539     1.6543      1.6545   5.533  < 2e-16 ***
urban                         0.2357     0.3837      0.3838   0.614 0.539145    
 
(conditional average) 
                            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)                  12.0109     0.6931      0.6931  17.328  < 2e-16 ***
closed_forest                 4.8906     0.8534      0.8535   5.730  < 2e-16 ***
cultivated                   12.1076     1.1800      1.1801  10.259  < 2e-16 ***
elevation_delta_scaled        0.5572     0.9383      0.9385   0.594 0.552689    
herbaceous_vegetation        11.4514     1.2644      1.2646   9.055  < 2e-16 ***
herbaceous_wetland           11.7865     3.6101      3.6106   3.264 0.001097 ** 
mean_elevation_scaled        -3.9688     1.5221      1.5223   2.607 0.009134 ** 
number_of_checklists_scaled  24.8715     0.7373      0.7374  33.729  < 2e-16 ***
open_forest                  17.4451     1.1501      1.1502  15.167  < 2e-16 ***
permanent_water              -0.9816     0.8237      0.8238   1.191 0.233477    
realmAfrotropic              -0.6229     1.4407      1.4409   0.432 0.665529    
realmAustralasia             -3.0481     2.1240      2.1243   1.435 0.151329    
realmIndomalayan             -4.2920     0.8641      0.8642   4.967    7e-07 ***
realmNearctic                -2.5242     0.9463      0.9464   2.667 0.007648 ** 
realmNeotropic               -3.9439     1.0244      1.0245   3.849 0.000118 ***
shrubs                        9.1539     1.6543      1.6545   5.533  < 2e-16 ***
urban                         0.5165     0.4215      0.4216   1.225 0.220493    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

blrHigherUrban_richness_sum <- model_summary('birdlife', 'richness', blrHigherUrban_richness_result)
blrHigherUrban_richness_sum

birdlifeHR_richness_plot <- blrHigherUrban_richness_sum[blrHigherUrban_richness_sum$model == 'full',]

birdlifeHR_richness_plot$explanatory <- factor(birdlifeHR_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
birdlifeHR_richness_plot$type <- 'Proportion of landcover\n5 km around site'
birdlifeHR_richness_plot$type[birdlifeHR_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
birdlifeHR_richness_plot$type[birdlifeHR_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
birdlifeHR_richness_plot$type[birdlifeHR_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'

birdlifeHR_plot = ggplot(birdlifeHR_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                 position=position_dodge(0.05)) +
  scale_y_discrete(
    limits = rev(levels(birdlifeHR_richness_plot$explanatory)), 
    labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
  theme_bw() +
  geom_vline(xintercept=0, linetype="dotted") +
  guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
  theme(legend.justification = "top")
birdlifeHR_plot
ggsave('percentage_of_regional_richness__output__birdlife_richness_higher_urban.jpg')

Saving 7.29 x 4.51 in image

Birdlife [Check pseudoreplication]: Difference in richness

locality_samples = read_csv('percentage_of_regional_richness__input__sampled_localities_min_distance_500m.csv')

── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
cols(
  LOCALITY_ID = col_character(),
  SAMPLE_ID = col_double()
)

table(locality_samples$SAMPLE_ID)

   0    1    2    3    4    5    6    7    8    9 
7391 7377 7388 7387 7395 7380 7382 7391 7382 7381 

plot_sample = function(sample_id) {
  birdlife_sample0 = birdlife[birdlife$locality_id %in% locality_samples$LOCALITY_ID[locality_samples$SAMPLE_ID == sample_id],]
  
  birdlife_sample0_richness_result <- model_average_richness(birdlife_sample0)
  birdlife_sample0_richness_summary <- model_summary('birdlife', 'richness', birdlife_sample0_richness_result)
  
  birdlife_sample0_richness_plot <- birdlife_sample0_richness_summary[birdlife_sample0_richness_summary$model == 'full',]
  
  birdlife_sample0_richness_plot$explanatory <- factor(birdlife_sample0_richness_plot$explanatory, levels = c('closed_forest', 'open_forest', 'shrubs', 'herbaceous_vegetation', 'herbaceous_wetland', 'permanent_water', 'cultivated', 'urban', 'mean_elevation_scaled', 'elevation_delta_scaled', 'realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic', 'number_of_checklists_scaled'))
  birdlife_sample0_richness_plot$type <- 'Proportion of landcover\n5 km around site'
  birdlife_sample0_richness_plot$type[birdlife_sample0_richness_plot$explanatory %in% c('mean_elevation_scaled', 'elevation_delta_scaled')] <- 'Scaled (0-1) average\nelevation 5 km around site'
  birdlife_sample0_richness_plot$type[birdlife_sample0_richness_plot$explanatory %in% c('realmAfrotropic', 'realmAustralasia', 'realmIndomalayan', 'realmNearctic', 'realmNeotropic')] <- 'Realm'
  birdlife_sample0_richness_plot$type[birdlife_sample0_richness_plot$explanatory %in% c('number_of_checklists_scaled')] <- 'Scaled (0 - 1) number\nof checklists'
  
  birdlifeSample0_plot = ggplot(birdlife_sample0_richness_plot, aes(y=explanatory, x=richness_birdlife_estimate, colour = type)) + 
    geom_line() +
    geom_point()+
    geom_errorbar(aes(xmin=richness_birdlife_estimate-richness_birdlife_error, xmax=richness_birdlife_estimate+richness_birdlife_error), width=.2,
                   position=position_dodge(0.05)) +
    scale_y_discrete(
      limits = rev(levels(birdlife_sample0_richness_plot$explanatory)), 
      labels = c('closed_forest' = 'Closed forest', 'open_forest' = 'Open forest', 'shrubs' = 'Shrubs', 'herbaceous_vegetation' = 'Herbaceous\nvegetation', 'herbaceous_wetland' = 'Herbaceous\nwetland', 'permanent_water' = 'Permanent\nwater', 'cultivated' = 'Cultivated', 'urban' = 'Urban', 'mean_elevation_scaled' = 'Mean elevation\nscaled', 'elevation_delta_scaled' = 'Elevation delta\nscaled', 'realmAfrotropic' = 'Afrotropical', 'realmAustralasia' = 'Austaliasian', 'realmIndomalayan' = 'Indomalayan', 'realmNearctic' = 'Nearctic', 'realmNeotropic' = 'Neotropical', 'number_of_checklists_scaled' = 'Sqrt number\nof checklists')) +
    theme_bw() +
    geom_vline(xintercept=0, linetype="dotted") +
    guides(colour=guide_legend(title="Predictor type")) + xlab('Increase in proportion of regional pool richness\n± Standard Error') + ylab('Predictor') +
    theme(legend.justification = "top")
  birdlifeSample0_plot
}

plot_sample('0')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample0.jpg')

Saving 7.29 x 4.51 in image

plot_sample('1')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample1.jpg')

Saving 7.29 x 4.51 in image

plot_sample('2')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample2.jpg')

Saving 7.29 x 4.51 in image

plot_sample('3')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample3.jpg')

Saving 7.29 x 4.51 in image

plot_sample('4')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample4.jpg')

Saving 7.29 x 4.51 in image

plot_sample('5')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample5.jpg')

Saving 7.29 x 4.51 in image

plot_sample('6')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample6.jpg')

Saving 7.29 x 4.51 in image

plot_sample('7')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample7.jpg')

Saving 7.29 x 4.51 in image

plot_sample('8')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample8.jpg')

Saving 7.29 x 4.51 in image

plot_sample('9')

Warning: comparing models fitted by REMLFixed term is "(Intercept)"

ggsave('percentage_of_regional_richness__output__birdlife_locality_sample9.jpg')

Saving 7.29 x 4.51 in image